/*!	 polynomial_root.h
**	 Polynomial Root Finder Header
**
**	Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
**	Copyright (c) 2007 Chris Moore
**
**	This package is free software; you can redistribute it and/or
**	modify it under the terms of the GNU General Public License as
**	published by the Free Software Foundation; either version 2 of
**	the License, or (at your option) any later version.
**
**	This package is distributed in the hope that it will be useful,
**	but WITHOUT ANY WARRANTY; without even the implied warranty of
**	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
**	General Public License for more details.
**
*/

#ifndef __SYNFIG_POLYNOMIAL_ROOT_H
#define __SYNFIG_POLYNOMIAL_ROOT_H

#include <complex>
#include <vector>

template < typename T = float, typename F = float >
class Polynomial : public std::vector<T> // a0 + a1x + a2x^2 + ... + anx^n
{
public:

    // Will maintain all lower constants
    void degree(unsigned int d, const T & def = (T)0)
    {
        resize(d + 1, def);
    }
    unsigned int degree()const
    {
        return this->size() - 1;
    }

    const Polynomial & operator+=(const Polynomial &p)
    {
        if (p.size() > this->size()) {
            resize(p.size(), (T)0);
        }

        for (int i = 0; i < p.size(); ++i) {
            (*this)[i] += p[i];
        }

        return *this;
    }

    const Polynomial & operator-=(const Polynomial &p)
    {
        if (p.size() > this->size()) {
            resize(p.size(), (T)0);
        }

        for (int i = 0; i < p.size(); ++i) {
            (*this)[i] -= p[i];
        }

        return *this;
    }

    const Polynomial & operator*=(const Polynomial &p)
    {
        if (p.size() < 1) {
            this->resize(0);
            return *this;
        }

        unsigned int i, j;
        std::vector<T> nc(*this);

        // in place for constant stuff
        for (i = 0; i < nc.size(); ++i) {
            (*this)[i] *= p[0];
        }

        if (p.size() < 2) {
            return *this;
        }

        this->resize(this->size() + p.degree());

        for (i = 0; i < nc.size(); ++i) {
            for (j = 1; j < p.size(); ++j) {
                nc[i + j] += nc[i] * p[j];
            }
        }

        return *this;
    }
};

class RootFinder
{
    std::vector< std::complex<float> >	workcoefs;

public:
    std::vector< std::complex<float> >	coefs; // the number of coefficients determines the degree of polynomial

    std::vector< std::complex<float> >	roots;

    void find_all_roots(bool polish);
};

#endif